SizeOverviewCombinedPlots

library(tidyverse)
library(here)
library(plotly)
library(htmlwidgets)

load(file = here("data/data_tidy.RData")) # load object data_tidy

BodyLengths <- tibble(Species = c("Ceriodaphnia sp.", "C. megalops", "D. cucullata", "D. curvirostris", "D. hyalina var. gellata", "D. hyalina var. lacustris", "D. hyalina", "D. longispina", "D. pulex", "D. magna", "Calanoid copepods", "Cyclops", "Bosmina coregoni", "B. longirostris", "Sida sp.", "S. crystallina", "Chydorus ovalis", "Eurycercus lamellatus", "Alona spp.", "A. quadrangularis", "Asplanchna", "Keratella spp.", "Ostracod", "Diaptomus"),
                      BodyLength = c(0.4636, 0.925, 0.8593, 1.58, 0.9064, 0.9064, 0.9064, 1.0347, 1.1656, 1.4214, 0.6860, 1.0369, 0.4230, 0.3637, 0.5075, 0.5075, 0.475, 1.975, 0.5235, 0.9679, 0.4747, 0.3495, 1.25, 0.9886)) # species and their body sizes

data <- left_join(data_tidy, BodyLengths, by = "Species") # add species body size to the dataframe

data <- mutate(data, SizeClass = case_when(BodyLength <= 0.6 ~ "Small (<= 0.6 mm)",
                                           BodyLength <= 1.0 ~ "Medium ( 0.6 < x <= 1.0 mm)",
                                           BodyLength > 1.0 ~ "Large (> 1.0 mm)"))
plotlist <- list()
lakes <- unique(data$LakeName)

for (i in 1:length(lakes)){
  Ldata <- filter(data, LakeName == lakes[i])
  
  if (nrow(Ldata) > 0) {
    
    summary <- Ldata %>% group_by(Species, SizeClass, Date) %>% summarize(sum_counts = sum(Counts, na.rm = TRUE),
                                                                          sum_total = sum(Total_individuals, na.rm = TRUE)) # Calculate counts per Date
    
    summary <- summary %>% group_by(SizeClass, Date, sum_total) %>% summarize(sum_counts = sum(sum_counts, na.rm = TRUE))
    
    summaryV <- Ldata %>% group_by(Species, Date) %>% summarize(sum_counts = sum(Counts, na.rm = TRUE),
                                                                sum_total = sum(Total_individuals, na.rm = TRUE))
    
    summaryV <- summaryV %>% mutate("Relative_abundance" = sum_counts / sum_total * 100) # Calculate relative abundance per Date
    
    summaryW <- left_join(summaryV, BodyLengths, by = "Species")
    
    WeightedMean <- summaryW %>% group_by(Date) %>% summarize(WeightedMean_mm = weighted.mean(BodyLength, Relative_abundance))
    
    #print(knitr::kable(WeightedMean, 
    #format = "markdown",
    #caption = paste0("The weighted mean for ", lakes[i], " in ", years[j])))
    
    summary <- summary %>% mutate("Relative_abundance" = sum_counts / sum_total * 100) # Calculate relative abundance per Date
    
    scaleFactor <- max(summary$Relative_abundance) / max(WeightedMean$WeightedMean_mm)
    
    result <- ggplot()+
    geom_col(data = summary, aes(x = Date, y = Relative_abundance, fill = SizeClass), color = "black", linewidth = 0.05)+
    geom_point(data = WeightedMean, aes(x = Date, y = WeightedMean_mm * scaleFactor), color = "darkmagenta")+
    scale_y_continuous(name = "Relative abundance (%)", sec.axis = sec_axis(~./scaleFactor, name = "Weighted mean in mm"))+
    theme_minimal()+
    theme(
        axis.title.y.right=element_text(color = "darkmagenta"),
        axis.text.y.right=element_text(color = "darkmagenta"))+
    labs(x = "Date", title = paste0("Relative abundance of each size class and the weighted means per date of ", lakes[[i]]))
    
    plotlist[[(length(plotlist)+1)]] <- result # add the plot to the list of plots
    
  } else {
    warning(paste0("No data for ", lakes[[i]]))
  }
}

htmltools::tagList(lapply(1:length(plotlist), function(x) { ggplotly(plotlist[[x]]) }))